# title: greene2005 estimation algorithm
# author: StableGenius
# date: 2020-01-22
# description: 使用maxLik暴力复现Greene2005提出的估计算法
# references:
## - GREENE W. Fixed and Random Effects in Stochastic Frontier Models[J]. Journal of Productivity Analysis, 2005, 23(1): 7–32.
## - GREENE W. Reconsidering heterogeneity in panel data estimators of the stochastic frontier model[J]. Journal of Econometrics, 2005, 126(2): 269–303.

# import packages
library(plm)
library(tidyverse)
library(frontier)
library(maxLik)

# sum grouped by individual
within_sum <- function(data, tdim){
  cumsum_t <- cumsum(tdim)
  cumsum_t <- c(0, cumsum_t)
  
  if (is.vector(data)){
    if (sum(tdim) != length(data)){
      stop("sum of tdim should be equal to length of vector!")
    }
    result <- numeric(length(tdim)) 
    for (index in seq_along(tdim)){
      result[index] <- sum(data[(cumsum_t[index]+1):cumsum_t[index+1]])
    }
  } else if (is.matrix(data)){
    if (sum(tdim) != nrow(data)){
      stop("sum of tdim should be equal to nrow of matrix!")
    }
    result <- matrix(NA, nrow = length(tdim), ncol = ncol(data))
    for (index in seq_along(tdim)){
      tmp <- data[(cumsum_t[index]+1):cumsum_t[index+1], ]
      if (is.vector(tmp)){
        result[index, ] <- tmp
      } else if(is.matrix(tmp)){
        result[index, ] <- colSums(tmp) 
      }
    }
  }
  
  result
}


# "true" fixed-effect model estimation
sfa_tfe <- function(formula, data, model = "product"){
  
  # check model specification: production or cost
  if (model == "product" | model == "production"){
    S <- 1
    ineffDecrease <- TRUE
  } else if (model == "cost"){
    S <- -1
    ineffDecrease <- FALSE
  } else {
    stop("Wrong model specification!")
  }
  
  # check panel data
  if (!"pdata.frame" %in% class(data)){
    stop("data set is not panel data!")
  }
  
  # dimension of panel data
  N <- pdim(data)$nT$n
  Ti <- pdim(data)$Tint$Ti
  
  # get variables
  y <- model.response(model.frame(formula, data))
  X <- model.matrix(formula, data)
  # drop intercept
  if ("(Intercept)" %in% colnames(X)){
    X <- X[, -1]
  }
  
  # parameters initialization
  fe_model <- plm(formula = formula, data = data, effect = "individual", model = "within")
  sfa_model <- sfa(formula, data, ineffDecrease = ineffDecrease)
  ## beta(no intercept)
  beta_0 <- coef(fe_model)
  K <- length(beta_0)
  ## sigma
  sigma_0 <- sqrt(coef(sfa_model)["sigmaSq"])
  ## lambda
  sfa_gamma <- coef(sfa_model)["gamma"]
  lambda_0 <- sqrt(sfa_gamma / (1-sfa_gamma))
  ## alpha parameters vector
  alpha_0 <- fixef(fe_model)
  # construct and rename
  parameters <- c(beta_0, sigma_0, lambda_0, alpha_0) 
  names(parameters) <- c(names(beta_0), "sigma", "lambda", paste0("id", 1:N))

  # log likelihood 
  logL <- function(param){
    param_beta <- param[1:K]
    param_sigma <- param[K+1]
    param_lambda <- param[K+2]
    param_alpha <- param[(K+3) : length(param)]
    
    if(param_sigma<0 | param_lambda<0)
      return(NA)
    
    ep <- (y - param_alpha - (X %*% param_beta)) %>%
      as.vector
    
    sum(log(2) - log(param_sigma) + log(dnorm(ep/param_sigma)) + 
          log(pnorm(-S*param_lambda*ep/param_sigma)))
  }
  
  # gradient of log likelihood
  logL_grad <- function(param){
    param_beta <- param[1:K]
    param_sigma <- param[K+1]
    param_lambda <- param[K+2]
    param_alpha <- param[(K+3) : length(param)]
    
    ep <- (y - param_alpha - (X %*% param_beta)) %>%
      as.vector
    
    p_star <- dnorm(ep*param_lambda/param_sigma)
    P_STAR <- pnorm(-S*ep*param_lambda/param_sigma)
    
    grad_beta <- colSums(X*ep)/param_sigma^2 + S*param_lambda/param_sigma*colSums(X*p_star/P_STAR)
    grad_sigma <- sum(-1/param_sigma + ep^2/param_sigma^3 + S*param_lambda/param_sigma^2*p_star/P_STAR*ep)
    grad_lambda <- -S/param_sigma * sum(p_star/P_STAR*ep)
    grad_alpha <- within_sum(ep/param_sigma^2 + S*param_lambda/param_sigma*p_star/P_STAR, Ti)
    
    c(grad_beta, grad_sigma, grad_lambda, grad_alpha)
  }
  
  # Hessian of log likehood
  logL_hess <- function(param){
    param_beta <- param[1:K]
    param_sigma <- param[K+1]
    param_lambda <- param[K+2]
    param_alpha <- param[(K+3) : length(param)]
    
    ep <- (y - param_alpha - (X %*% param_beta)) %>%
      as.vector
    
    p_star <- dnorm(ep*param_lambda/param_sigma)
    P_STAR <- pnorm(-S*ep*param_lambda/param_sigma)
    pp <- p_star/P_STAR^2
    block <- -P_STAR + param_lambda^2/param_sigma^2*P_STAR*ep^2-S*param_lambda/param_sigma*p_star*ep
    pp_block <- pp*block
    colsum_pp_block <- colSums(X*pp_block)
    
    hess_beta <- -t(X) %*% X/param_sigma^2 + (S * param_lambda^2/param_sigma^2) * (t(X) %*% diag(pp*(param_lambda/param_sigma*P_STAR*ep-p_star*S)) %*% X)
    hess_beta_sigma <- -2*colSums(X*ep)/param_sigma^3 + S*param_lambda/param_sigma^2*colsum_pp_block
    hess_beta_lambda <- -S/param_sigma * colsum_pp_block
    hess_beta_alpha <- (-X/param_sigma^2 + S * param_lambda^2/param_sigma^2 * X * pp * (param_lambda/param_sigma*P_STAR*ep-p_star*S)) %>%
      within_sum(Ti)
    hess_sigma <- sum(1/param_sigma^2 -3*ep^2/param_sigma^4 - 2*S*param_lambda/param_sigma^3*p_star/P_STAR*ep+
                        S*param_lambda^2/param_sigma^4*pp*(param_lambda/param_sigma*P_STAR*ep^3-S*p_star*ep^2))
    hess_sigma_lambda <- -S/param_sigma^2*sum(pp_block*ep)
    hess_sigma_alpha <- (-2*ep/param_sigma^3 + S*param_lambda/param_sigma^2 * pp_block) %>%
      within_sum(Ti)
    hess_lambda <- S/param_sigma^2*sum(pp*ep^2*(param_lambda/param_sigma*P_STAR*ep-S*p_star))
    hess_lambda_alpha <- (-S/param_sigma*pp_block) %>%
      within_sum(Ti)
    hess_alpha <- (-1/param_sigma^2 + S*param_lambda^2/param_sigma^2*pp*(param_lambda/param_sigma*P_STAR*ep-p_star*S)) %>%
      within_sum(Ti)
    
    result <- matrix(NA, nrow = K+2+N, ncol = K+2+N)
    result[1:K, 1:K] <- hess_beta
    result[K+1, K+1] <- hess_sigma
    result[K+2, K+2] <- hess_lambda
    result[(K+3):(K+2+N), (K+3):(K+2+N)] <- diag(hess_alpha)
    
    result[K+1, 1:K] <- hess_beta_sigma
    result[1:K, K+1] <- hess_beta_sigma
    
    result[K+2, 1:K] <- hess_beta_lambda
    result[1:K, K+2] <- hess_beta_lambda
    
    result[K+2, K+1] <- hess_sigma_lambda
    result[K+1, K+2] <- hess_sigma_lambda
    
    result[1:K, (K+3):(K+N+2)] <- t(hess_beta_alpha)
    result[(K+3):(K+N+2), 1:K] <- hess_beta_alpha
    
    result[K+1, (K+3):(K+N+2)] <- hess_sigma_alpha
    result[(K+3):(K+N+2), K+1] <- hess_sigma_alpha
    
    result[K+2, (K+3):(K+N+2)] <- hess_lambda_alpha
    result[(K+3):(K+N+2), K+2] <- hess_lambda_alpha
    
    result
  }
  
  maxLik(logLik = logL, 
         grad = logL_grad, 
         hess = logL_hess, 
         start = parameters, 
         method = "BFGS")
}